Landau level states on a topological insulator thin film 
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We analyze the four-dimensional Hamiltonian proposed to describe the band structure of the 
single-Dirac-cone family of topological insulators in the presence of a uniform perpendicular mag- 
netic field. Surface Landau level(LL) states appear, decoupled from the bulk levels and following 
the quantized energy dispersion of a purely two-dimensional surface Dirac Hamiltonian. A small 
hybridization gap splits the degeneracy of the central n = LL with dependence on the film thick- 
ness and the field strength that can be obtained analytically. Explicit calculation of the spin and 
charge densities show that surface LL states are localized within approximately one quintuple layer 
from the surface termination. Some new surface-bound LLs are shown to exist at a higher Landau 
level index. 
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I. INTRODUCTION 

Insulating materials with topologically protected sur- 
face states known as topological insulators (TIs) arc a 
matter of great current interest 1-3 . The surface metal- 
lic states in this new class of materials is characterized 
by Dirac-like quasiparticle dispersion, and a one-to-one 
correspondence between momentum and spin quantum 
numbers of the single-particle states thus representing 
an extreme form of spin-orbit coupling. Both these as- 
pects have been confirmed for the first time in Bi^Sbi-a; 
family 4 of topological insulators by ARPES 5 and STM 6 
studies. 

More recently, a lot of experimental efforts has been 
given to the synthesis and characterization of Bi2Se3, 
Bi2Te3, and Sb2Te3 7-14 following the prediction of their 
topological behavior 15,16 , due to their simple surface 
band structure consisting of a single-cone Dirac spectrum 
centered at the T-point and a relatively large band gap. 
Topological insulators of the singlc-Dirac-conc family in 
the thin-film form has been synthesized by a number of 
groups 11-13 . Theoretically, the thin-film TIs bear close 
analogy to another heavily studied topological material, 
i.e. graphene 17 . For instance, the well-known pair of 
valley-degenerate Dirac bands of graphene becomes the 
top and bottom surface Dirac bands of TIs with finite 
thickness. Perpendicular magnetic field quantizes the 
surface Landau levels (LLs) with the energies that scale 
with the LL index n as ±-y/n in TIs as well as in graphene. 

Previous treatments of the magnetic field effect on 
TI surface started from the two-dimensional (2D) Dirac 
Hamiltonian focusing only on the surface electronic states 
and ignoring the bulk states altogether 18,19 . These meth- 
ods relied on first projecting the bulk Hamiltonian to the 
surface, obtaining the 2D Dirac model, then including the 
field effect by way of Peierls substitution. In another vein, 
several recent papers theoretically examined the proper- 
ties of a thin slab of TI in which the bulk and surface 
electronic states are treated on an equal footing 20-22 in 



the absence of the magnetic field. It is thus natural to 
consider how the magnetic field effect plays out for a 
thin film geometry of TI, following the spirit of solving 
the bulk Hamiltonian adopted in Refs. 20-22. In fact, 
an attempt of precisely this sort has been made in a re- 
cent paper by Liu et al. 16 Here, the authors solved the 
4x4 tight-binding Hamiltonian with the Peierls sub- 
stitution for the magnetic field and even including the 
Zeeman field coupling. Wc point out in this paper that 
the method adopted in Rcf. 16 does not treat the sur- 
face and bulk electronic states simultaneously, and as a 
result the bands arising from surface LLs penetrate into 
the bulk LL states, while physically such overlapping of 
energy levels will not occur. 

Our approach follows closely the spirit of zero- field case 
studied in Refs. 20-22 and takes care of the boundary 
conditions properly. Some parts of our report are tech- 
nical, dealing with the characteristic equation resulting 
from the boundary conditions and the methods of solv- 
ing them. Several physically meaningful results follow 
from our analysis. First, hybridization of the zeroth-LL 
states localized on the top and the bottom surfaces for a 
sufficiently thin sample is shown to manifest itself as the 
splitting of the degeneracy of zeroth-LL states with the 
gap magnitude that can be calculated analytically. The 
zeroth-LL gap size oscillates with the film thickness. Our 
finding naturally extrapolates a similar observation of the 
gap oscillation observed previously 20 to finite magnetic 
field. Interestingly, we find that a new kind of surface- 
bound LL states appear for highcr-LL indices where the 
conventional surface LL band of ~ \fn variety has merged 
into the bulk continuum. Justification of the new sur- 
face LLs is made on the basis of careful numerical study 
and an approximate analytic solution of the characteris- 
tic equation. Properties of the bulk single-particle states 
for higher-LL indices are examined in detail. Finally, 
both charge and spin density profiles of the surface LLs 
at low LL indices along the thickness of the sample are 
explicitly worked out. 



In Sec. II we formulate the LL problem based on the 
4x4 Hamiltonian proposed previously for Bi 2 Se 3 -family 
of topological insulators. Boundary conditions are im- 
posed on the two surface layers for a thin-film geometry 
and characteristic equations are derived in Sec. III. In 
Sec. IV several physical results are shown and its rele- 
vance to recent STM are discussed. Summary of results 
and an outlook is given in Sec. V. Technical discus- 
sion for the new surface-bound LLs can be found in the 
Appendix. 
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II. FORMULATION 

The 3D tight-binding Hamiltonian proposed as a mini- 
mal model for singlc-Dirac-cone family of TIs first in Rcf. 
15 and detailed in Ref. 16 is 



FIG. 1: (color online) (a) Surface energy spectra without mag- 
netic field when Ek = (red) and et 7^ (blue) . (b) Bulk and 
surface energy dispersions in the absence of magnetic field and 

£1- = 0. L z /l z — 3000 was used. 
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trices r are introduced and p± = p x ± ip y are momentum 
operators. The upper and lower indices in the basis set 
refer to the parity and spin quantum numbers for the p z 
orbitals of Bi or Se atoms, respectively. It was shown 15 
that e(p) and M(p) depend on the momentum p as 



e(p) = C + D iP 2 z + D 2 (pI+pI), 
M(p) - M -B x pI-B 2 (pI+pI] 



(2) 



Values of the various constants can be found in Refs. 
15,16,20-22. In our paper all the material parameters 
are re-scaled in terms of the one mass scale Mq. Two 
length parameters emerge as a result, l z = Ai/M and 
l± = A2/M0, each characterizing the length scale within 
the plane and perpendicular to it. With the material 
parameters given in Rcf. 15 they read lj_ = 14.64A and 
l z = 7.9 A. We use them as the measure of length in each 
direction. All equations can be cast in dimensionless form 
as well as the two functions e(p) and M(p) which now 
become (following the parameterization of Ref. 15) 



fc+u = ( Et x + k z + iMt v \ v, 
k—V= ( Et x — k z + iMr y J u, 



M 
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(4) 



with E as the energy, k± = k x ± ik y as the momentum, 
and a z and a± are two material constants. They read 
a, = 0.58 and aj_ = 0.94 in the paramctrization of Rcf. 
15. 

As our interest lies in the case of finite thickness L z 
for the z-direction the above equation will be deformed 
as ik z — > A z 20 ~ 22 . Due to the boundary conditions at 
z = ±L z /2, surface state solutions appear 20-22 . Figure 
1(a) shows the difference in the surface energy spectra 
when the diagonal energy eu is turned on/off. As stressed 
earlier we will suppress the diagonal energy e^ and work 
with the particle-hole symmetric model in the following 
section where we consider the magnetic field effect. Fig- 
ure 1 (b) shows the surface state energy together with the 
bulk energy as a function of the transverse momentum 
k±. The edge state energy dispersion is precisely linear 
in \k±\ for large |fcj_| but opens an exponentially small 
hybridization gap at k±_ = 0. The gap at the T-point is 
given by 22 



e(p) = -0.024 + 0.075p 2 + 0.3265p^ 
M(p) = 1 - 0.58p 2 - 0.94pi. 



(3) 



Cocfficient-by-coefficient, expressions in e(p) are smaller 
than the ones in M(p). In this study, we will ignore 
e(p) for calculational simplicity and restore particle-hole 
symmetry of the spectrum as a consequence. 

The four-dimensional single-particle eigenstates can be 
constructed in terms of two, two-dimensional spinors u 
and v. For an infinite medium one can write the eigen- 



state as ij) 



r X, where x = 



is a 4-component 



constant spinor to be determined by solving 



8a 



where a, /3 are 



sin(/3L,)|, 



1 , VS 

a = - — , p = 



2a 2 



2a 2 



(5) 



(6) 



This result will be generalized in the following section 
to the nonzero magnetic field, with the revised meaning 
for the gap as the energy difference of symmetric and 
anti-symmetric combinations of zcroth-Landau levels lo- 
calized to top and bottom surface layers. 



III. LANDAU LEVELS 



Magnetic field H = (0, 0, H) perpendicular to the 
slab modifies the momentum operator p > p + A = 
(p x — Hy,p y ,p z ) in the Hamiltonian. A pair of canonical 
operators 



are introduced such that [A, A^] = 1. The magnetic 
length (measured in units of l±) appears as Ir = 1/VH, 
as well as the guiding center j/q = lnk x . Relation to the 
physical held strength i? p h ys in Tesla is 



H = lleH phys /h ~ 3.25 x Kr 3 # phys /[T]. 



(8) 



Taking k + = -{V2/l H )A^ and fc_ = -(y/2/l H )A, the 
eigenvalue equation for a slab with perpendicular mag- 
netic field becomes 



I 

-j= [et x + i\ z + iM^ u, (9) 



A f u = — y= \Et x - iX z + iMfjTyj v, 

Av -- 



with several new definitions (N = A^A) 



a H 



^i, M fi = l + a z \ 2 z - aH (N+\y (10) 



'H 



The rest of this section is concerned with the solution of 
this equation, together with the boundary conditions at 
the two terminations z = ±L z /2. 

The structure of the equation invites for a solution of 

is the n-th Landau level (LL) oscillator wave function 
centered at y = y$. By substituting the ansatz to Eq. 
(9) we get 23 



the form u 



a, 



and v 



where 
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V2 
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1 + a z X z — anin 



1 



(11) 



We can parameterize the spinor solution u and v sat- 
isfying Eq. (11) in the following form 



tan 9 = 



E + [i 



, tany> 



M n _i+ M 



2n/l H 



(13) 



Here (i means 



lia H = Ml - E 2 - A 2 + a H M n 



2n 

1 a 



(14) 



The eigenvalues are fixed up by the relation fi 2 = E 2 + A 2 
which reads when /i is explicitly written out 



Mi-E A -\{ 



a H M n 



2n 

' H 



a 2 H (E 2 + X 2 z ). (15) 



This is the desired characteristic equation for the energy 
E. 

Being eighth-power in X z , one can find eight different 
A 2 's for a given energy. We call them aXb as in the non- 
magnetic case 20-22 , with a = ± and b = 1,2,3,4. There 
are thus eight independent solutions of the same energy 
E for a given LL index n and the guiding center yo, 



Xnab(y-yo) 



-la sin Ub cos ip^fpn-i 

COS 9b COS ipb4>n-l 

cos 9 b sin (fib4>n 
ia sin 9b sin (pb<fin 



(16) 



Taking the linear combination among the eight states 
gives out the most general eigenstate before the boundary 
condition is imposed as 



rl>nk. (x, y~y Q ,z) = e ik * x ^ A ab e aXbZ X nab{y - W)).(17) 

ah 

To facilitate the further solution, the coefficients A a b 
can be classified into symmetric (A a b = Ab) and anti- 
symmetric (A a b = aAb) types. For the symmetric case 
the boundary conditions at z = ±L z /2, tp n k x (x,y — 

2/o, W 2 ) = i>nk x { x >y - 2/0, --W 2 ) = °> can be satisfied 
if we require that Ab obey 



>J A(,sinh(A(,L z /2)sin6';,sin( / 3( ) = 0, 

b 

2_] Ab sinh(AhL z /2) sin 9b cos ipb = 0, 

b 

>J Ab cosh(\bL z /2) cos 9b sin tpi, = 0, 
b 

\] Ab cosh(XbL z /2) cos 9b cos tpj, = 0. (18) 



9„_1 COS If 



-i sin9 
cos 9 



UTn^) 



with the two complex angles ((9, ip) fixed by 



A nontrivial solution exists provided the characteristic 
equation of the above 4x4 matrix is zero. With the aid 
of Eq. (13) this condition can be expressed as 



\aL z 



' Ai A 2 tanh ^-^ tanh ^p A 3 A 4 tanh ±f± tanh 

- {E n +^n,l)(E n +fJ-n,2) {E n +fJ-n,3)(En+fJ-nA) 

■ Ax A 4 tanh ^* tanh ^- A 2 A 3 tanh ^^ tanh ^^ 

(E n + jJ. n ,%){E n + Mn,3) 



(.E„ + [l n ,l)(E„ + Mn,4j 

Ai A 3 tanh ^^ tanh ^^ 

(£„ +/i„,l)( £: «+Mr l ,3) 



A2A4 tanh 



A 2 £* 



tanh 



A 4 L 3 



(-E n + /U n ,2)(-E'rj+Mn,4) 



{fin,l—fin,2+a z (Xl - A2))(/U n ,3-Mn,4+Q1 Z (A3 - A4)) 
J(/Vl-Mn,4+a z (Ai _ ^l))(A i n,2-A i r l ,3+a 2 (A2 - A3)) 



(A*n,i-Mn,3+a z (A? - A3))(^„, 2 -A i n,4+a z (A2 - A4)), 



(19) 



where M n . b = 1 + a z Aj - a H (n + 1/2), and [i n ^a H = M% b - E* - \\ + a H M n<b + 2n/l 2 H . 



The case of anti-symmetric coefficients A a b = aAb 
can be handled by interchanging smh(XbL z /2) and 
cosh(A&L z /2) in Eq. (18), and replacing tanh(Ab-L z /2) 
by coth(A(,i z /2) in Eq. (19). Equation (19) and its 
anti-symmetric counterpart can be solved numerically for 
given n, giving out simultaneously surface and bulk en- 
ergy solutions in the presence of the field H . When L z 
becomes large both tanh(A&L 2 /2) and coth(AbL z /2) tend 
to the same value and we will have a pair of degenerate 
states for each energy, each state being localized either 
at the top or the bottom surface and not coupled to the 
opposite layer. 

The zeroth-LL n = requires a separate treatment. In 
this case u is identically zero, and v T = (co,do) is found 
from solving 



^(Er.-iX.+iM^) (2 



0, (20) 



with Mq = 1 + a z \ 2 z — an/2. For a given energy Eq, 
Eq = Mq — X\ results in four different A z 's, a\ with 
a = ± and b = 1,2. Similar to Eq. (16) one can assume 
the spinor solution for n = 



^ 

cos 6» 6 
\ia sin 9b 



(21) 



where Ob is given by 



tan #& 



A h 



E + Mo.b' 



(22) 



and 



M , b = l + a z X 2 b -a H /2, 

\ b = — = — (l-2a z + a' ± a z 
V2a z v 

-(-1) Vl - 4«, + 2a H a z + AE 2 a 2 z ) * .(23) 

A linear combination 



i>ok x (y -yo) = ^2 A abe aXbZ xoab(y - s/o) 



(24) 



can be formed with the boundary conditions at z = 
±L z /2. Again assuming symmetric (A a b = Ab) and anti- 
symmetric (A a b — aAb) coefficients separately and denot- 
ing the corresponding energies by Eq and Eq, we have 



and 



A2 E +Mq ] i 
JlEfi + Moa 



X 2 E Q A +M 0A 



tanh^Jj^ 



tanh 



tanh 



A 2 L 3 



A 2 i = 



XiE*- 



-M , 



tanh^i 



(25) 



(26) 



One can easily show from Eqs. (25) and (26) that Eq = 
-E s 

Cjq . 

This completes the derivation of the full energy spec- 
tra and eigenstates for a thin-slab geometry of TI model 
with the perpendicular magnetic field. In the following 
section we discuss several physical results obtained from 
the analysis of the solution. 



IV. PHYSICAL RESULTS 

Figure 2 shows the dependence of surface and bulk en- 
ergies on the LL index n, for a sufficiently large thickness 
L z . The numerical results remain consistently similar for 
L z larger than about ten times l z . A most surprising as- 
pect of the numerical analysis is the existence of three 
distinct branches of surface-localized states, labeled as 
(I), (II), and (III) in Fig. 2. 



(i) 



The behavior of the first surface branch E„ 
ably close to the formula: 



E { n s) ~ ±V2^/l H = ±V2nH. 



is remark- 



(27) 



This is exactly what is expected of the purely two- 
dimensional Dirac Hamiltonian with the Fermi velocity 



v F = I (equal to M l±/h = A 2 /h w 6.2 x 10 5 m/s in 
physical units) 15 . Restoring all physical units, the sur- 
face LLs occur at 



E {s) = ±- 



-h 



Ih„ 



'2n. 



(28) 



Physical magnetic field H p h ys = 11T results in the mag- 
netic length lH phys = \Z^/eH p i iys ~ 100A and the energy 
levels ±58y/n mV. Indeed the spacing in the n = and 
n = l LL peaks were found to be about 40 ~ 50 mV 13 . 

Using the physical magnetic field H p \ iys = 10T we find 
that at n > n c ps 12 the first branch of surface LL begins 
to merge with the bulk spectrum (Fig. 2). Here n c corre- 
sponds to the Landau level index for which the surface LL 
begins to touch the bottom of the bulk band. A sharper 
criterion to determine n c can be drawn by keeping track 
of the eigenvalues Af, for the surface-bound LLs. With 
increasing n, one of the four Ab's forming the surface LL 
cigenstate has its real part decrease and eventually touch 
zero at n = n c . This signals the mixture of an extended 
state in the wave function just as the surface LL merges 
with the bulk continuum. Recently, the number of sur- 
face LLs that can be resolved in the tunneling spectra of 
STM 13 was shown to be about 12, consistent with our 
estimate of n c . For n > n c , the surface branch no longer 
exists independently of the bulk LL, but rather seems to 
form the bottom of the bulk band as depicted in Fig. 2. 




FIG. 2: (color online) Landau level energies for L z /l z =3x 
30 3 and H p hy S =10T, showing both surface(sky blue square) 
and bulk states(black square). Three surface branches are 
labeled (I) through (III) with analytic fits shown as red solid 
curves to y2n/lu (branch I) and y/A ± 8 — Bn (branches II 
and III). A and B values are derived in the Appendix and 
a small offset 5 is used to fit the numerical results. When 
Hp hys =10T, A = 0.918, 5 = 0.068, B = 0.037, respectively. 

A recent paper by Liu et al. also computed the bulk 
and surface LLs based on the model Hamiltonian for 
Bi2Sc3. In their Fig. 7 it appeared as though the surface 



LLs can exist well inside the bulk spectra as an indepen- 
dent branch. We believe this is an artifact of their cal- 
culation not taking care of the boundary conditions pre- 
cisely. Once the boundary conditions at z = ±L z /2 are 
handled properly, the correct energy profile for n > n c is 
the one in which the surface-localized wave functions are 
hybridized with the extended states to form a "hybrid" 
state. To confirm this assertion, we have made a careful 
analysis of all the A;, values for eigenstates with ener- 
gies both at the bottom of, and deep inside the bulk for 
7i > 7i c . While the details are too tedious to report here, 
we can say with certainty that states forming the bulk 
LL are typically a linear combination of solutions with 
real Af, (localized to surface) and some with purely imag- 
inary Xb (extended). See Eq. (17) for a general definition 
of the cigenstate. Only for the three surface branches (I) 
through (III) is it possible to get all Ab's of the eigenstate 
being real and the wave function completely localized. 

The existence of extra two surface branches, labeled 
(II) and (III) in Fig. 2, is unexpected. They begin to 
appear at n w 4 and n ~ 8 respectively for H p h ys = 10T. 
We have confirmed their existence for L z /l z as small as 
10 and as large as 3000. Due to the insensitivity of their 
features to surface thickness, we can first of all conclude 
that the extra surface modes are bound to one particular 
surface and not hybridized with the other one. To further 
confirm that these branches are genuine, we have carried 
out an approximate analytic treatment valid at large LL 
index n and infinite thickness L z and indeed found that 
two extra branches exist. Details of this analysis are 
given in the Appendix. 
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FIG. 3: (color online) Hybridization gap energies for H p h ys — 
0T, 10T, and 50T with varying thickness L z . 

Hybridization effect mixes the two degenerate n = 
LLs previously associated with each surface layer and 
opens a gap. We have derived the n = surface LL ener- 
gies analytically for symmetric (Eq ) and anti-symmetric 
(Eq ) combinations as 



E* 



and Eq 



Aa(H) 



(l-a ± H)sm[P(H)L z }e 



-a(H)L z 



,(29) 



-Eq 3 . The gap is defined as A = |£^ 



i-0 4 



2\Eq\. For practically available field strengths 
where H < 1, a{H) and /3(i/) in Eq. (29) are 



a(H) 



1 
2o~ 



/3(F) 



V4a 2 -l-4a 2 a ± F 
2o~ ' 



(30) 



They reduce exactly to a and /3 coefficients obtained in 
Eq. (6) as H — > 0. The gap still exhibits an oscillatory 
decay similar to the gap at the T-point without magnetic 
field. In Fig. 3 we compare the energy gaps for zero-field 
and for -ff p h ys = 10T and 50T. The similarity of their In- 
dependence is a strong clue that the origins of the gaps 
are the same. Ignoring the small field-induced shift, the 
gap can be 



A « 7M e-W[9.2A] w 2e- L >^- 2A ^V. 



(31) 



It gives a value w 10 meV for a seven quintuple-layer 
thin film and may well be resolved as two split n = LLs 
in a careful STM spectroscopy study. Currently avail- 
able thin-film STM study was done on 50 quintuple-layer 
sample 13 . In Rcf. 21 it was argued that the oscillation 
in the sign of the hybridization gap under zero magnetic 
field marks the transition between topologically trivial 
and non-trivial insulator phases. If this is so, our cal- 
culation seems to reveal that well-defined Dirac-like LLs 
exists regardless of the thickness and the sign of the gap, 
implying that changes in the topological character of the 
thin-film TI will not be revealed by examination of the 
surface LLs alone. 
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FIG. 4: (color online) nth-LL wave function density p n ,c(z) 
and the spin density p„ >s (z)ior L z = 60A. 

The charge (c) and spin (s) densities of each n-th sur- 
face LL wave function can be defined as 



Pn,s(c){ z ) — 



J dxdy iptr s{c) 4> r , 



J dxdydz ijjnipn 
T s ■--- diag(l, 1,-1,-1), 
r c = diag(l,I,l,l). 



(32) 



Figure 4 shows results for a few surface LLs with small LL 
index n. All surface LLs arc localized to within one L of 



the termination, or within about one quintuple layer. As 
one can see from Fig. 4(b), the zeroth-LL is completely 
spin-polarized, / po, s {z)dz = — 1, while other higher sur- 
face LLs are nearly spin-quenched, f p n >o,s(z)dz w 0. 
The zeroth-LL has only the lower two elements of the 
four-component spinor \ take nonzero values, which re- 
fer to the amplitudes for Bi and Se states of spin- 1 (See 
text following Eq. 1). There are two n = LL in the so- 
lution, and both of them are fully spin-4.-polarized. The 
origin of the spin polarization is the analogue of the sub- 
lattice polarization of the n — LL in graphene 17 . The 
difference is that the two valley n = Landau levels 
occupy the opposite sublattices, so that the overall sub- 
lattice symmetry is restored. 

Here, by contrast, both top and bottom surface LLs 
give the same spin polarization. The reason is that for 
the top surface Dirac states the magnetic field is pointing 
out of the bulk but the bottom surface states experience 
the field pointing into the bulk, so that effectively the 
sense of the field direction is also reversed between the 
two surface layers. By reversing the field direction from 
+z to — z one will generate n = of spin-'f polarization. 
As a result a thin slab of TI subject to quantizing mag- 
netic field creates two n — LLs which are completely 
spin polarized. Such spin-polarized surface layers are de- 
tectable by Faraday or Kerr rotation experiments 15 . 



CONCLUSION 



We showed how to derive the Landau level solution 
for a slab geometry of the topological insulator based on 
the four-band model 15,16 . Previoius approaches were to 
first project the zero-field bulk Hamiltonian to the sur- 
face, then using the Peierls substitution to address the 
magnetic field effect 16 ' 19 . Our strategy by contrast is to 
introduce the Peierls substitution directly into the bulk 
Hamiltonian and use the boundary conditions appropri- 
ate for a slab geometry. The obtained surface Landau 
level energies are in good accord with those obtained from 
the surface Dirac Hamiltonian, and we conclude that sur- 
face projection and the Peierls substitution can be imple- 
mented in any order with the same physical spectrum. 

A dramatic departure of the present Dirac LL problem 
with an analogous one posed by the graphene system 1 is 
that the surface LLs are eventually bounded by the bulk 
spectra, and one has to face the issue what will happen 
to the surface LLs as they begin to merge with the bulk 
continuum. We addressed such a question numerically 
and analytically in this paper, with a prediction for the 
existence of new surface-bound LLs appearing at higher- 
LL indices. Detection of the predicted new surface modes 
presents an interesting challenge for the future surface- 
sensitive measurements on TI materials. 
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Here m b is a constant, to be determined later. This gives 
for M n , 



M.„ 
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Appendix A: Analysis of New Surface Modes 

After some trial and error, we find that the following 
ansatz describe the numerically found surface modes (2) 
and (3) with good accuracy. 



We can also make an ansatz for the surface energy mode 
of the form 



Ei =A-Br 



(A3) 
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with two undetermined positive coefficients A and B. In- 
(Al) serting Eqs. (Al) through (A3) into Eq. (15) gives 
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The terms in • • • have subleading order in n than the ones 
shown. Assuming a sufficiently large n we require that 
the two sides of the equation cancel out at each order in 
n. From the equality of n 2 , n, and -^/n-order terms we 
obtain the three following equations: 
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Upon solving them we obtain 
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Further consideration of sub-leading corrections finally 
yield a splitting of A into two branches responsible for 
(II) and (III) in Fig. 2. Rather than going into the 
complicated sub-leading order analysis, we can simply 
split A into two branches by writing A±8 with d chosen to 
fit the two branches in Fig. 2 while A itself is completely 
determined from the parameters such as a z and a±. It 
is shown that both branches (II) and (III) match quite 
well the ansatz for energy, Eq. (A3). 

We can also discuss the stability of the new surface 
branches by recalling the numerical value a z = 0.58, 
and Oh = 2a±/l H = 1.84/l H . It follows that pos- 
itive (a H m b ) 2 is possible if 2 - (0.58)" 1 - 1M/1 2 H = 
0.276 - 1.84/1% > 0, or if 1% > 6.66. Returning to physi- 
cal length scales, this implies the magnetic length greater 
than ~38A, or the magnetic field strength less than 46T. 
We then expect that surface modes (II) and (III) should 
co-exist with the more familiar mode (I) inside the bulk 
gap for typical laboratory magnetic field ranges. 
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